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Abstract. These notes summarize basic concepts underlying numer- 
ical relativity and in particular the numerical modeling of black hole 
dynamics as a source of gravitational waves. Main topics are the 3+1 
decomposition of general relativity, the concept of a well-posed initial 
value problem, the construction of initial data for general relativity, 
trapped surfaces and gravitational waves. Also, a brief summary is 
given of recent progress regarding the numerical evolution of black 
hole binary systems. 



1 Introduction 

The theory of general relativity (GR) has enchanted generations of physicists through 
its conceptual and mathematical beauty, which is expressed in the simplicity and 
geometric nature of the Einstein equations, 

Gab = &TrnT ab , 

relating the curvature of spacetime to its matter content. It is in particular the 
connection of physical concepts such as the equivalence principle to a geometric 
description of nature, which is stressed in introductions to the field. General rela- 
tivity makes many exciting physical predictions, such as the big bang, black holes, 
or gravitational waves. The theory also provides a seemingly inexhaustible supply 
of mathematical challenges, and many deep mathematical insights have been gained 
in trying to understand the physical content of the Einstein equations, such as the 
positive mass theorem [1] or the nonlinear stability of Minkowski space [2] . The cur- 
rently emerging fields of gravitational wave astronomy, or, more generally, general 
relativistic astrophysics, are however changing this picture on a rapid timescale. A 
large international effort is underway to establish a network of gravitational wave 
detectors, such as LIGO |3l4j . GEO [516] and VIRGO [7j, some of which are already 
taking data at design sensitivity, and first publications have set new upper limits 
on the radiation from several sources (see e.g. |8|9|10j ). 

The field where the contact between mainstream general relativity and new 
applications is perhaps most intense, is numerical relativity (NR). Here one tries to 
systematically explore the solution space of the Einstein equations with numerical 
methods. One of the most advertised goals is to model the inspiral and collision of 
two black holes, and the resulting gravitational wave signals, in order to provide 
templates for gravitational wave data analysis. Many other applications are no less 
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exciting, e.g. the study of cosmological singularities {i.e. the predictions of general 
relativity about the very early universe), or critical collapse, where one probes the 
regime of extremely high curvatures created from the collapse of matter fields, i.e. a 
regime where classical general relativity will ultimately lose its validity and quantum 
effects will play a role. 

While one of the most timely challenges for numerical relativity is the establish- 
ment of a "data analysis pipeline" , connecting analytical calculations of the early 
inspiral phase with numerical simulations and gravitational wave searches in actual 
detector data, on a technical level such simulations lead to challenging questions not 
only regarding the physics of general relativity but also in applied mathematics, 
high performance computing and software engineering. 

In order to set up reliable and efficient algorithms that predict detector signals 
from "first principles" -type numerical solutions of the Einstein equations it is essen- 
tial that the individual steps are clearly understood, and subtleties or open issues 
be made transparent. Already the concept of gravitational wave emission from some 
source is highly nontrivial in general relativity, and touches upon the global struc- 
ture of spacetimes. Setting up the Einstein equations as an initial value problem 
and solving them numerically raises several new problems in the theory of partial 
differential equations and their numerical solution. The complexity of the Einstein 
equations requires a sophisticated approach toward writing software in order to 
produce reliable results. Efficiency is paramount, when one actually aims to make 
the connection to data analysis, where large parameter studies are expected to be 
required. Understanding optimal choices in this regard requires a careful analysis of 
the interplay of the continuum equations, numerical methods and their implemen- 
tation on current generations of high performance computing equipment. 

A key idea when one wants to study the solutions of the Einstein equations in a 
systematic way, is to write the Einstein equations in the form of a well posed initial 
value problem. That this is possible, is not trivial, and an important test of the the- 
ory. Like the other fundamental theories of interactions (the electroweak interaction 
and quantum chromo-dynamics), general relativity can be interpreted as a gauge 
theory (in this case of the diffcomorphism group), and initial value formulations 
thus contain constraint equations. These constraint equations, which are associated 
with the diffeomorphism invariance, and the general subtleties of diffeomorphism 
invariance, like the absence of a fixed spacetime background which is essential in 
most areas of physics, create challenging problems for the mathematical treatment 
of the classical theory, its quantization, and for any treatment with approximation 
methods, such as numerical simulations. In a sense, the root of the problems to 
quantize gravity, is the same as for many problems to solve the Einstein equations 
numerically. 

Discretizing GR can be approached in very different ways: 

— One approach, that also has found applications in quantum gravity, is the direct 
discretization of the geometry (e.g. in Regge calculus [TTj). 

— A related idea is to use "discrete differential forms" , in analogy of a very suc- 
cessful method in electromagnetism, see e.g. [12j . 

— The mainstream approach, which I will follow in these notes, is to view the 
Einstein equations as a standard PDE (partial differential equations) problem. 

All of the approaches however share the fundamental problem that due to the 
presence of the rather complicated constraints we have no direct access to physical 
degrees of freedom. As a consequence approximate solutions of constrained systems 
may not manifestly satisfy the constraints, it turns out that exponential drift off 
the constraint surface is typical, as is illustrated by the Maxwell example of eq. 
(J9j> . It is possible to trade evolution equations for constraints when constructing the 
solution - a minimal number of (hyperbolic) evolution equations would correspond 
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to the local degrees of freedom (i.e. two for vacuum gravity). This may improve 
the quality of numerical evolutions, but there is no guarantee: since there are more 
equations than evolution variables, one still needs to worry whether the evolution 
equations that are not used to construct the solution will be (manifestly) satisfied. 

The good news for numerical relativity is of course that in contrast to pertur- 
bation techniques numerical simulations do not rely on the smallness of physically 
relevant parameters, and can thus be applied in very general situations. Carefully 
crafted numerical simulations provide error estimates and convergence results. A 
particular challenge for the field, for instance, is to clarify the range of validity of 
post-Newtonian (PN) expansions [13], and significant progress has been made in 
this direction |14|15|16|17|18j . 

2 Gravitational Waves 

Like electromagnetism, general relativity predicts waves: aspherically accelerated 
masses emit gravitational radiation. Both Maxwell's and Einstein's theories have 
two local degrees of freedom, which correspond to the two polarization states of 
waves in the theory. While in electromagnetism waves carry no monopole moment, 
in relativity there are no monopoles and no dipoles, to leading order wave emission 
is thus quadrupolar. 

Under "everyday conditions" , including the conditions in a gravitational wave 
detector, the amplitude of gravitational waves is extremely small, and it is sufficient 
to consider the linearized theory. Consider the approximation g ao = r) a b + h a b for 
weak gravitational fields, where r\ a b is the metric of Minkowski spacetime, and h a b 
is a small perturbation. Using the definition 

h a b = h ao — —rjabh c c 

one can show that 

Uh ab = 

if the gauge condition 

dh ab /dx b = 

is satisfied. The tensor h ao thus satisfies the standard wave equation on a flat back- 
ground. It can be shown that a further gauge condition can be adopted to also set 
the trace of h a b to zero, n ab h ao = 0. In this " transverse-traceless (TT) gauge", 
the metric perturbation h ab corresponding to a wave traveling in the ^-direction 
of a Cartesian coordinate system can be written as a linear combination of two 
polarizations h + and h x : 

h i:j = h + (e + )ij + h x (e x )ij, (1) 

where e +jX are the basis tensors 

(e + )y- = iiij - y t yj , and (e x )ij = x l y J + Xjjji , (2) 

and a; and y are the unit vectors in the x and y directions respectively. For an 
extensive introduction to the linearized theory of gravitational wave, including the 
principles of their detection, see [TH] ■ 

When modeling sources of gravitational waves in full general relativity it turns 
out that a notion of gravitational waves becomes highly nontrivial: because of dif- 
feomorphism invariance there is no background spacetime available, on which grav- 
itational waves can be defined in an unambiguous way. Such a definition is only 
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possible asymptotically for isolated systems. When the geometry "flattens out" to 
approach Minkowski spacetime at large distance from the sources, then Minkowski 
spacetime can be used as a suitable background, see sec. O 

The weak field approximation can also be used to compute the radiation from 
sources, e.g. accelerated bodies. For velocities » « c and for wavelengths A much 
larger than the size of the system, a weak field calculation yields the loss of energy 
of a system of massive bodies as 

dE 

dt 5c 5 ^ \ dt 3 J 1 ' 

with Qij — f g(xiXj — ^5ijr 2 )d 3 x is the mass quadrupole moment. Eq. §5§ is known 
as the quadrupole formula (for the energy loss of the system). The radiation power 
scales with the sixth power of the frequency of the system. Due to the weakness of 
gravity, expressed in the factor , thus only systems of astrophysical dimensions 
- large masses moving at a significant fraction of the speed of light - generate 
significant amounts of gravitational radiation. While the motion of the earth around 
the sun only generates around 200 Watt of gravitational radiation, a close neutron 
star binary (100 km separation, 100 Hz) would generate approximately 10 45 Watt. 
In the inspiral and collision of two black holes, approximately 4 % of the total energy 
is radiated away in the form of gravitational waves - for the collision of two stellar 
size black holes with ten solar masses each, this would already be around 0.8 solar 
masses! Indirect confirmation of the predictions of general relativity for the energy 
loss of a binary system due to the emission of gravitational wave has been possible 
by measuring the tightening of the orbit of the Hulsc- Taylor binary pulsar [20|21j . 

Two essential properties of gravitational waves for astrophysical observations are 
that they are not shielded by interstellar masses, and that detectors are sensitive 
to amplitude rather than intensity. The first means that not only will gravitational 
wave observations render possible the direct observation of phenomena that have 
not been observable so far, but we will also be able to see regions that have been 
hidden behind dust clouds. Gravitational waves represent the bulk motion of fast 
compact objects instead of incoherent superposition of many atoms, which is typi- 
cal of electromagnetic observations. The observation of coherent wave trains makes 
gravitational wave observations similar to detecting sound waves, rather than to 
producing images, and in particular results in a direct sensitivity to amplitude. 
Consequently the sensitivity scales with 1/ distance instead of 1/ distance 2 as for 
electromagnetic observations. This makes it much easier to observe objects at ex- 
treme distances, and also means that any improvement in detector sensitivity by 
some factor A yields a factor A 3 in event rates. 

3 The Initial Value Problem for General Relativity 
3.1 The Maxwell Equations 

GR is a classical relativistic field theory much like Maxwell's theory of electromag- 
netism, but it is nonlinear, and its dynamics involves the very structure of spacetime 
itself. We first discuss the structure of the Maxwell equations as an example. The 
Maxwell equations can be elegantly and geometrically formulated in a 4-dimensional 
way as 

dF = 0, d*F = 4tt * J, (4) 

where F is the Faraday tensor of electromagnetic field strengths, and J is the 
4-dimensional current. For practical purposes it is useful to solve the Maxwell equa- 
tions as an initial value problem, that is to specify suitable data at some instant 
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of time, and then evolve them to the future. The formulation as an initial value 
problem allows a systematic approach to constructing solutions: First one needs to 
classify permissible initial data, then the future and/or past development of such 
data is determined uniquely by evolution equations. The Maxwell equations ^ are 
not written in this form, in particular it is not immediately obvious what permis- 
sible initial data would be, and whether time evolution of such data prescribes a 
unique solution. 

A standard way to proceed is to prescribe a foliation of spacetime, e.g. by speci- 
fying a time function i(x^), such that the surfaces St, defined by t = const., form a 
1-parameter foliation of spacelike surfaces. One can then perform a geometric split 
of 4-dimensional tensorial quantities into 3-dimensional objects (i.e. objects defined 
in the tangent and co-tangent space of St), and consider the equations for these 
3-dimensional objects. It will prove useful to define a timelike unit normal to the 
surfaces S t , 

Tin ■= • 

V c iV c i 

The timelike unit normal can be used to define projection operators onto timelike 
(vertical) and spacelike (horizontal) directions. It is easy to check that n a rib projects 
onto the direction of the timelike unit normal, and the induced metric 

Kb = 9ab + n a n b (5) 

projects onto the tangent space of the spacelike hypersurface St, where g a b denotes 
the (Lorentzian) spacetime metric, and h ab then is the induced Riemannian (positive 
definite) metric on S t . A covariant derivative D a on S t that is compatible with the 
metric h ab , D a h bc = 0, can be defined as 

D rpai...a r — h c 'h ai h ar h bl h b ' B X7 , r T a t— a 'r 

Excellent pedagogical discussions of this geometric splitting procedure can be found 
in [22123] , 

For the Maxwell equations one can now define 

E a = F ab n», B e = ^F cd h c a h\e abef n f , 

the electric and magnetic field as projections of the field strength tensor F an , as well 
as the charge density p = J a n a and current j b — J a h ab . On Minkowski spacetime, 
and also using the standard slices of Minkowski spacetime, the Maxwell equations 
separate into the following two groups of equations in terms of the electric and 
magnetic fields: 

D a E a =^p, D a B a = 0, (6) 

and 

d t E a = e abc d b B c ~ 47T Ja , (7) 
d t B a = -e abc d b E c . (8) 

Since Eqs. (J6j) do not contain time derivatives, they have to be interpreted as con- 
straint equations that restrict the space of allowed initial data, while (17|8j) are evo- 
lution equations. It is not hard to show that if the constraints are satisfied initially, 
they will be satisfied for all times, furthermore, the evolution equations determine 
the time evolution in a unique way. Without these two facts, the initial value de- 
scription of the Maxwell equations, eqs. ((6][11) would make little sense. Because of 
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their linearity, the discussion of the initial value problem, and the solution of the 
constraints are tremendously simpler for the Maxwell equations than for the Ein- 
stein equations, but the essential structure is the same. Some remarks are in order: 
Counting degrees of freedom in eqs. (JB][H]), we find that there are 2 local degrees of 
freedom (6 first order evolution equations minus 2 constraint equations makes 4 first 
order in time or 2 second order in time equations), these correspond to the two po- 
larization states of the electromagnetic field, ft is clear that while for the continuum 
solution the constraints propagate (they are satisfied at all times by virtue of the 
evolution equations if they are satisfied at the initial time), this is not necessarily 
true for a numerical implementation. Discretization error may give rise to a growth 
in the constraints, which may lead to instability or are least spoil the physical valid- 
ity of the approximate solution. We will see that the same problem appears for the 
Einstein equations. In electromagnetism, many approaches have been developed to 
deal with this issue, see |24j for a comparison. The most standard approach uses 
a formulation of the Maxwell equations in differential forms language, and directly 
translates the differential forms concepts to the discrete level [25]. This way, the 
constraints can be solved manifestly in the discrete problem. 

The problems of preserving the constraints of the Maxwell theory become more 
pronounced in curved spacetime, and already in curved slices of Minkowski space- 
time, where the constraint propagation equations become 

CnD.E 1 = -KDiE\ C n D t B l = —KD i E i (9) 

where C n is the Lie derivative along the timelike unit normal of the hypersurface, 
and K = D a n a is the trace of the extrinsic curvature of the hypersurface (see below). 
Here the sign of K is chosen such that K > corresponds to expansion, and K < 
to collapse. Clearly, in the collapsing case, initial violations of the constraints will 
be amplified. While a simple redefinition of the fields (densitizing them) solves the 
problem in the Maxwell case |26j , this still provides a valuable toy model for similar 
problems in numerical relativity [27] . 

3.2 The Einstein Equations 

The Einstein equations G a t = SnuTab can be written as an initial value problem 
in formal analogy with the previous discussion of the Maxwell equations. There 
is much freedom in how to do this, and the route we will take is to very briefly 
introduce the most traditional way, and then to discuss its issues with regard to 
numerical relativity. 

First, we need to make a careful choice of the topology of our spacetime manifold 
M, since the Einstein equations allow solutions that violate causality in a global 
way, i.e. geometries where there exists no spatial hypersurface St that uniquely 
determines the geometry at all earlier and later times. A necessary requirement for 
causality [23] is that the topology of the spacetime be £ x I, where £ is a three- 
dimensional manifold, and / is a (possibly infinite) interval. We can then choose 
coordinates such that the line element takes the form 

ds 2 = ~a 2 dt 2 + lab {dx a + /3 a dt)(dx b + (3 b dt). (10) 

In particular we have thereby singled out a time function t. 

We can now write the Einstein equations in terms of 3-dimensional objects: the 
induced metric h a t (the definition in (JTDJ) coincides with the definition ([5]) , the shift 
vector f3 a and the lapse function a. The coordinate basis vector (d/dt) a can be 
written in terms of the unit normal n a , and the lapse and shift as 




= an a + [3 a , n a f3 a =Q. 
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A reduction to first order in time of the field equations can be obtained by intro- 
ducing the extrinsic curvature (or second fundamental form) 

K ab := -Cnjab = — ('jab - D a f3b - D b /3 a ) , (11) 
2 2a 

C n is again the Lie derivative with respect to the vector field n a . Making a particular 
choice of adding constraints to the evolution equations one arrives at a "standard 
form" of the 3+1-decomposed Einstein equations, which has dominated numerical 
relativity for several decades, and is commonly referred to as the ADM-equations, 
and called the York- ADM equations here, since these equations can be viewed as a 
variant [28] of an evolution system discussed by Arnowitt, Deser and Misner in [29] : 

D b K\ - D a K\ = 0, 3 R+(K a a ) 2 ~K ab K ab = 0, (12) 

K ab = -D a D b a + (3 c D c K ab + K cb D a p c - K bc D a f3 c + a( 3 R ab + K c c K ab ), (13) 

where for simplicity matter terms have been set to zero. The first two equations 
can be interpreted as constraints analogous to the Maxwell divergence constraints, 
while the third equation together with the definition of the extrinsic curvature 
(fTTjl forms a first order in time evolution system. The Bianchi identity V a G ab = 
implies that the constraints propagate, that is they are satisfied by virtue of the 
evolution equations at all times if they are satisfied initially. For a much more 
detailed discussion of the 3+1 decomposition of the Einstein equations see [30] . 

The big advantage of the ADM equations is their relative simplicity. The problem 
for numerical relativity, which has only been realized after several decades, is that 
the "free evolution problem" , where the constraints are only solved initially and 
only the evolution equations are used to construct the solution at later time, is only 
weakly hyperbolic in the language of hyperbolic partial differential equations, and 
the initial value problem therefore ill-posed (with the exception of certain subclasses 
of initial data such as spherical symmetry, which has added further to the confusion) . 
The issue is discussed in detail in [31] . 



3.3 BSSN: the workhorse formulation for the binary black hole problem 



For a long time, the standard choice of variables for writing the Einstein equa- 
tions was based on the York- ADM equations 32J. It is known now, however, that 
the hyperbolic subsystem thus obtained is not well posed; specifically, it leads to a 
weakly hyperbolic set of equations (see e.g. [31]). A formidable industry of creating 
improved evolution systems has produced numerous alternatives, the most popular 
being the "BSSN family" 33 34 35 . This formulation is characterized by introduc- 
ing a contracted connection term as a new variable, a conformal decomposition of 
the metric and extrinsic curvature variables, and adding constraints to the evolution 
equations. 

Detailed discussions of well-posedness for the BSSN family have been given 
by Gundlach and Martin-Garcia [36|37|38] , The set of evolved variables are the 
conformally rescaled unimodular three- metric 7^ , the logarithm of the conformal 
factor </>, the trace of the extrinsic curvature K, the conformally rescaled traceless 
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extrinsic curvature Ay, and the contracted Christoffcl symbols r z : 

= -^log(det 7 y), (14) 

%=e- A % j , (15) 
K = 1 ij K ij , (16) 

Ay = e- 4 HK tJ - pijK), (17) 

r = q k f k . (is) 

As is usual, we will adopt the convention that indices of densitized quantities (de- 
noted with a tilde) are raised and lowered with the conformally rescaled three-metric 
7y . The introduction of new variables leads to new constraints, one differential and 
two algebraic: 

G = f l - ^ k r l jk =0, S = det 7 y - 1 = 0, A = Al = 0, (19) 

which are again propagated by the evolution equations. 

The standard Hamiltonian and momentum constraints of general relativity take 
the form [39] 

H = e -4 * (R - 8D j D^ - &{jy<l>)(Dj4>) \ + \k 2 - A VJ A 11 - \jJC, 

Mi = QAU{D 3 4>) - 2A(D i( f>) - ^(DiK) + ^{D 3 A kl ). 

The BSSN evolution equations, which are obtained from the Einstein equations 
by using the definitions (|14p and making a standard choice for adding constraints, 
become 



r 



aK 
~~6~' 

C n K = -fl'fl,a + oAyA^' + 
£„Ay = -e-^(D l D J a) TF + a (V^(i?y) TF + KA l3 - 2A lk A k J 
C n r l = -2{d 3 a)A^ + 2n{l -, A^ - ^{d 3 K) + 6A«(c>^)), 

where A is covariant derivative associated with 7 y, C n — dt — Cp is the Lie deriva- 
tive along the unit normal, 7y denotes the trace-free part of a tensor Ty. The 
Ricci curvature in terms of the BSSN variables takes the form 

-0./2)^ k aidk% + i k(l d 3) r k + r k r [lj)k + 2j lm f k (l f 3)km + i lm r k m r kl3 . 

The algebraic constraints are typically solved at every time step, e.g. by setting 
7« -» Tydet^- 1 / 3 and Ay Ay - iA iro7 ^r' ro . 

Currently, the standard choice for evolving the lapse function for BSSN-cvolutions 
is given by the so-called Bona-Masso family of slicing conditions: 

dtot = -aKf(a), (20) 
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in particular the choices / = 1, corresponding to harmonic time slicing, and / = 
2/ a, which is usually termed "1+log" slicing, see [40] . 

In order to obtain long-term stable numerical simulations, it is equally important 
to construct a suitable shift vector field (5 % . Here we report on evolutions where 
we evolve the shift vector according to a Gamma- freezing prescription [40]. A key 
feature of this particular choice is to drive the dynamics of the variable F % towards 
a stationary state. As a "side effect" this choice creates a coordinate motion that 
drags the black holes along an inspiral orbit. A crucial effect of this method is that 
the resulting coordinate motion which corresponds to the naive physical intuition 
reduces artificial distortions in the geometry, which otherwise could easily trigger 
instabilities 

In summary, the solution procedure for the equations is as follows. First, we 
specify free data motivated by quasi-equilibrium arguments, then solve the 9 com- 
ponents of the constraint equations (TL, M.i, Q % , A, S) to obtain initial data for the 
17 evolution variables (</>, jij, K,Aij,r l ). The evolution system is completed by 
specifying evolution equations for the four gauge quantities (a,(3 l ), which yields a 
hyperbolic system that is second order in space and first order in time, and which 
determines the evolution of all 21 components of the "state vector" describing the 
geometry of spacetime. 

In finite difference codes for the solution of nonlinear hyperbolic equations, it 
is common practice to add artificial dissipation terms to all right-hand-sides of the 
time evolution equations, schematically written as 

Otu^dtu + Qu. (21) 

Such terms damp out high-frequency noise, e.g. as produced by mesh-refinement 
boundaries, and can be necessary to guarantee numerical stability for nonlinear 
problems [41] . The typical form of this terms is the Kreiss-Oliger dissipation oper- 
ator [41] (Q) of order 2r 

Q = a(-h) 2r - 1 (D + ) r P (D^y/2 2r , (22) 

for a 2r — 2 accurate scheme, with a a parameter regulating the strength of the 
dissipation, and p a weight function that we typically set to unity. 

Discretization in space is performed with standard second-, fourth- 42 43] or 
sixth-order [44] accurate stencils. In particular, symmetric stencils are used, with 
the exception of the advection terms associated with the shift vector, asymmetric 
upwind stencils are used. Time integration is performed by standard Runge-Kutta 
type methods, in particular 3rd and 4th order Runge-Kutta and second order ac- 
curate three-step iterative Crank-Nicholson integrators as described in [31], where 
Courant limits and stability properties are discussed for the types of equations used 
here. 



4 Partial differential equations: Well-posedness and numerical 
stability for initial value problems 

An excellent and seemingly trivial starting point for a discussion of numerical ap- 
proximations is the model problem F(x, y) = 0. An important issue in the context 
of approximations is the sensitivity in the dependence of a solution y on an input pa- 
rameter x. It is useful to define a condition number, which quantifies the worst possi- 
ble effect on y when x is perturbed. Consider the perturbed eq. F(x+Sx, y+Sy) = 0, 
and define 

K = SUP \MUM\ 
5* I Ml/I Ml ' 
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If K is small we call the problem well conditioned, if K is large ill conditioned, and 
if K (y) = oo ill-posed or unstable. As a starting point for numerical relativity we 
clearly need the continuum problem to be well-conditioned - in GR this is far from 
trivial, e.g. because this will rely on a judicious choice of coordinate gauge. 

We call an evolution problem well posed if a unique solution exists (in a gauge 
theory such as general relativity this requires a gauge choice) and depends continu- 
ously on the initial data. The latter condition is usually expressed as the condition 

|K*)ll<^*|K0)||, (23 ) 

where the exponential term ensures robustness with respect to lower order terms, 
and the constants a and K can be chosen independently of the initial data. Note 
that this condition is only required local in time, since global in time solutions 
may not exist in a nonlinear theory - singularities may form! An obvious crucial 
question is which norm ||.|| is appropriate for defining well-posedness of a certain 
type of differential equation. It turns out that for first order in space and time 
systems the standard L2-norm is sufficient. Since the Einstein equations take the 
form of second differential order equations in the metric, a complete reduction to 
first order may seem artifical. Note that there is an important difference between 
first order reductions in space and in time: Reduction to first order in time leads to 
new evolution equations, but reduction to first order in space leads to new evolution 
and constraint equations. This enlargement of the phase space may not only reduce 
computational efficiency, but can also give rise to instabilities or pathologies not 
inherent in the original problem. 

First order in time formulations provide a "normal form" for ordinary differential 
equations and are thus also convenient for systems of PDEs, which in numerical 
analysis are often discussed in terms of the "method of lines". In this approach 
first only space is discretized, and time left continuous. PDEs are thus converted to 
coupled systems of ordinary differential equations. From a physical point of view, 
first order in time formulations are attractive, e.g. because they are most easily 
integrated into a Hamiltonian formulation. 

The concept of well-posedness translates straightforwardly to the concept of 
numerical stability for discrete iterative problems. Consider a simple stable iterative 
algorithm 

v n+1 = Q(t n ,v n )v n : \\v n \\ < Ke^Wv^W Vw°, 

where v the solution vector and Q a matrix. The stability criterion allows e atn 
growth, but excludes e an growth, i.e. for differential equations exponential growth is 
allowed in the continuum problem, but resolution dependent growth of the numerical 
algorithm is excluded. 

To illustrate the analysis of numerical stability for a simple ODE problem, con- 
sider the standard textbook example 

y' = Xy, y(0) = y a 

and solve numerically with the forward Euler algorithm, 

Vn+i = Vn + hy' n = y n + h\y n . 

The stability criterion yields |y„ + i|/|2/„| = |1 + h\\, and the algorithm is unstable 
for h > — 2/A. For A > the continuum solution grows exponentially and even 
stable algorithms will suffer from ill-conditioning. 

4.1 The wave equation as a toy model for hyperbolic equations 

The wave equation provides a standard example to illustrate well-posedness for 
evolution equations and to introduce different notions of the associated technical 
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concept of hyperbolicity. Starting with the second order form in ID, <pj t — c 2 <fi jXX , 
we can obtain a mixed first/second order version, 

<j>,t = ctt, 7T it = C(j} <xx 
or a complete first order reduction 

<j>,t = C7T, (j) tX =1p, K,t = Clp, x , ijjj = CTT.x, 

where = "0 now plays the role of a constraint which is preserved by the evolution 
equations: 

d t (<j>,x ~ ip) = d x d t 4> - d t i> = d x 7T - d x ir = 0. 

Note that the evolution equation for <f> decouples, and we may focus on the system 
of equations for 7r and which has the form 



dtu = Ad x u, u = {it, if}}, A 



c 
c 



We call the matrix A, or generally the coefficients of the highest spatial derivatives, 
the principal part of the system of partial differential equations. The characteristic 
variables v± — ip ± tt correspond to the eigenvectors of A and therefore satisfy 
advection equations, u± — ±cu iX , where the eigenvalues of A are ±c. We can easily 
construct the general solution in Fourier space: 

11^ = TTo cos(ci u>) + iipo sin(ci w), ip u = ipo cos(ctuj) + i 77o sm(ctui). 

The norm ||w|| 2 = J |77| 2 + l^l 2 can easily be checked to satisfy ||w(£)|| = ||u(0)||, 
well-posedness can thus directly be read off from the general solution, and numerical 
stability can be discussed in an analogous way. General theorems can be used to 
show stability against perturbations with lower order terms |45) . 
The system with principal part 

A, =00' (24) 

which we also will refer to as the 1- Jordan Block model, shows very different features: 
A' has real eigenvalues - (1,1), but does not have a complete set of eigenvectors 
(characteristics) and can not be diagonalized. The solution is U u = (u,v), u = 
Lot$muj(t + x) , v — sinw(i + x). Consider data with u(0) = and frequency to - 
in terms of the L2 norm we now get 



ll^(o)|| 



= Vl + t 2 uj 2 , (25) 



which does not satisfy our criterion for well-posedness. It turns out, that alternative 
norms can be chosen which render this problem well posed, but generic lower order 
perturbations convert the frequency dependent linear growth to frequency depen- 
dent exponential growth [45j . The choice of an appropriate norm is thus a subtle 
problem, and robustness with respect to lower order terms a crucial criterion. 



4.2 First order systems 



Consider a linear first order system with constant coefficients: 

u = Ad x u + Bu + C, 
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where u is interpreted as a multi-component object. It is easy to see that for a single 
equation, the function A corresponds to the propagation speed of the wave. For a 
system of equations, the eigenvalues of the matrix A clearly can again be interpreted 
as propagation speeds if the matrix A is diagonalizable. The eigenvectors are then 
called the characteristic variables. General theorems allow to restrict analysis of 
well-posedness and stability to B = C = [45]. The formal solution in Fourier 
space is 



In order to study whether the solutions can be bounded by initial data and thus 
whether the system of equations allows a well-posed initial value problem, we need 
to evaluate e IujA . Intuitively, no problems will arise if A is diagonalizable with real 
eigenvalues, the solution will then be purely oscillatory. 

Well-posedness and stability can indeed be discussed in terms of eigenvalues 
(characteristic speeds) and eigenvectors (characteristic variables) of A, (again we 
refer to [45] for an excellent textbook presentation): 

— If all eigenvalues (speeds) arc real, the system is called (weakly) hyperbolic, and 
is well posed in absence of lower order terms in an appropriate norm. 

— If a complete set of eigenvectors exists (the characteristic variables span the 
solution space), the system is called strongly hyperbolic, and admits a well posed 
initial value problem. 

— If the system is strongly hyperbolic and admits a conserved energy it is called 
symmetrizable (symmetric) hyperbolic (strongly hyperbolic implies symmetriz- 
able in ID). 

As an example consider the York-ADM equations (|11I12I13[) in ID with lapse 
function a — \/det g (a standard choice), linearized around Minkowski space: 



Computing the Jordan normal form of the matrix A characterizing the principal 
part of the first order system yields 



All characteristic speeds are real, but there are 2 Jordan blocks, the eigenvectors of 
the characteristic matrix thus do not span the solution space and the system is only 
weakly hyperbolic. Nevertheless, many physics results have been obtained with the 
York-ADM system in spherical symmetry, and excellent test results are obtained 
in the ID test suites of [46j . This is explained by noting that decoupling the xx 
components from the yy and zz components leads to well-posed systems and good 
numerical results. This happens for the "gauge wave" metric eq. (|26j) . a linearized 
wave on flat background: 



u = IujAu 



u = 





J(A) 



/00 000\ 

-1 

-11 

-100 

10 

11 

\0 1/ 



ds 2 = -dt 2 + dx 2 + (1 + b)dy 2 + (1 - b)dz 2 , b = Asm 




but also in spherical symmetry, where again the yy and zz components (or, say 
and if ip) of the metric are not independent. 
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4.3 Second order in space, first order in time systems 

First order in time, second order in space systems are typically used in astrophys- 
ically oriented numerical relativity codes, because of their simplicity and a vague 
expectation of better accuracy than first order reductions, which has at least par- 
tially been confirmed in [31] , Such systems are also attractive because they arise 
rather naturally in a Hamiltonian context. The issue of well-posedness for such 
systems has however only been clarified very recently 47 48 49 50 38 31j, for much 
of the history of NR. it was widely believed that a clean route to well-posedness 
requires the standard theory of first order hyperbolic systems, see e.g. |45j for an 
excellent textbook reference. 

Again, the wave equation provides an excellent simple example of what is going 
on. The first order in time, second order in space wave equation is 

d t (f>(t, x) = II(t, x), d t II(t, x) = d xx (/){t, x). 

Consider the family of solutions [51] 

4>{x,t) — sin(wx) cos(wi), 77(x,<) = — lo sin(wx) sin(urf) 

with initial data 4>o(x) = sin(wx), IIq(x) = 0. Varying u> in the initial data, the L2 

norm of the solution at time t, J ^ (\4>\ 2 + \n\ 2 ) dx, can be made arbitrarily large 
with respect to the initial data (whose norm is independent of lo), this contradicts 
well-posedness in L,2- Since numerical codes implementing these ill-posed equations, 
as well as the York- ADM evolution equations (fT3|) which had been argued to lead 
to an ill-posed initial value problem before, did not show blowup, it was suggested 
|51) that well-posedness may not be crucial for numerical relativity codes. 

The issue is resolved as follows: since all the solutions for the wave equation are 
actually bounded and show no growth, ill-posedness in L2 is not a sign of pathology 
in this case - the problem rather is the inappropriate use of the L2-norm. Translating 
the first order norm that one gets from introducing 

X = d x <b, 

which leads to a well posed, symmetric hyperbolic problem, to the second order 
variables, and using the norm implied by this translation 

/ 2 "(|77| 2 + |a^| 2 )dx 
Jo 

does in fact yield well-posedness, as can be shown by pseudo-differential reduction 
(i.e. solving explicitly in Fourier-space) [47j . Note that on the discrete level an 
ambiguity regarding the discrete norm needs to be resolved: it is not clear how the 
derivative should be discretized, or whether this matters at all. This ambiguity is 
resolved in |31j , see below. The reason that a code that implements the first order in 
time, second order in space wave equation does not show any pathological behavior 
is thus simple: the system is well posed, and as shown in |31j can be discretized 
stably in a straightforward way. The York- ADM equations are however indeed only 
weakly hyperbolic and therefore lead to an ill-posed initial value problem. However, 
the growth is only linear, and is easily dismissed without careful convergence tests, 
or when artificial dissipation is added. Using the norm that would correspond to a 
first order-reduction of the system notions of well-posedness and hypcrbolicity can 
then indeed be defined for mixed order hyperbolic systems [4714 8 49 50 3Tj. 
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A simple "normal form" for mixed order hyperbolic systems has been introduced 
in [31]: 

d t u = Pu u-f U ) P -f A ^ + B . C . 

Evolved variables are split into two types: u are differentiated twice and v are not. 
Not all second order in space systems can be written in this form (a simple example 
is the heat equation, u t = u xx , which is parabolic). Now consider the second order 
principal symbol P' and a matrix E, which is the principal part of the associated 
first order system: 



oo o 

o A n c 

D nn G r ' 



The main result on the continuum level of the equations is that if E is diagonalizable 
in a regular way, the initial value problem is well posed: 

Kt, Oil <Ke at \\u(Q r )\l H| 2 ^ [\u\ 2 +J2\dM 2 + \v\ 2 d d x, 

J »=i 

where the norm is again the straightforward translation of the L 2 norm of the 
associated first order system. This result is equivalent to results that had previously 
been obtained by other authors in 47 49 37J. The result may seem trivial, but 
has ended a long controversy in numerical relativity about whether the evolution 
systems typically used in mainstream numerical relativity are well posed, or even 
whether their well-posedness can be discussed in a mathematically rigorous way. 
Recent results along the lines sketched here prove well-posedness results for the 
evolution systems and coordinate gauge conditions that are actually used in large 
scale binary black hole evolutions [52], which had seemed far out of reach only a 
few years ago. 

The appropriate norm for the discrete case is 

d 

U^\\Al, D+ ^\Hl + h\\l + Y.\\ D ^ u \\l- 

i=l 

The ambiguity of discrctizing the derivative in the norm is resolved by using the 
first order forward (-D+) or backward (-D-) spatial derivative in combination with 
standard centered second or fourth order differencing. Using centered derivatives 
in the norm is shown not to be robust with regard to perturbations in lower order 
terms. 

Note that for the continuum version of the matrix E the frequency lo can be 
factored out, which is not generally the case for the discrete version of E, this would 
only hold in general if the second derivative is really also the derivative of the first 
derivative on the discrete level. For some evolution systems in general relativity 
certain standard discretizations are thus problematic, as discussed in [31j . 



4.4 Stable and well-posed is not enough 



Clearly, well-posedness and numerical stability are not sufficient to guarantee suc- 
cessful numerical simulations, since exponential growth or blowup in finite time, 
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which is allowed for well posed problems, is typically not tolerable for numerics 
(unless the timescale of the growth is sufficiently smaller than the timescale of in- 
terest in the solution). A typical trick to get rid of exponential growth is a change 
of variables, but this is in general difficult for tensors. As an example, it is a stan- 
dard result of textbooks on ordinary differential equations that sufficiently regular 
ODE initial value problems are always well-posed, solutions may however blow up 
in finite time: 

y' = ^y 2 ,y(Q) = yo -> y(x)=yo/(xy -i). 

In numerical relativity, exponential growth of the continuum solution already 
appears in seemingly trivial test problems. Consider a "gauge wave" metric - the 
flat metric of Minkowski space in coordinates that correspond to a traveling wave 
of coordinate distortion. This metric provides a challenging test for most evolution 
codes and is part of a numerical relativity test suite [46] : 

ds 2 = -Hdt 2 + Hdx 2 + dy 2 + dz 2 , H = H(x -t) = l-A sin (^^j- —J ■ (26) 

This line element represents Minkowski spacetime with a nontrivial choice of har- 
monic time slicing (V a V Q t = 0). One analytic source of rapidly growing error is 
the instability of flat space on T 3 . Another problem is the existence of a family of 
harmonic, exponential gauge modes corresponding to H — > e xt H , with arbitrary A. 
This problem can be modeled by a nonlinear wave propagating on Minkowski space 

1331 : 

^"dadb'P - ^(da^d^) = = $ V ab d a d b \0g$. 

Imposing periodic boundary conditions, one evolves on a 3-torus, i.e. there are no 
boundaries. For <P > the Cauchy problem is well-posed. In addition to the solution 
<P = 1 + F(t — z), F > —1, the system also admits the solutions 

$x = e xt {\ + F{t~z)) 

for arbitrary A. Numerical errors will excite the growing modes, which eventually 
dominate the signal. As pointed out in [53], an obvious solution is to use \og<P as 
an evolution variable, but the example just models a situation that does arise in 
numerical relativity, where it is not clear how to take the logarithm of the metric. 

For a problem that is actually ill-posed, the situation will be much worse than for 
our gauge-wave toy problem: typically higher frequencies will correspond to faster 
growth (i.e. larger a,K in the estimate ([23])). Since better resolution of the grid 
does in particular allow higher frequencies, improving the resolution will in general 
lead to a numerical solution that grows faster. This is analogous to an unstable 
numerical scheme - in such a case, the discretized equations do not lead to a well 
posed problem. 



5 Energy, momentum, and radiation 

For astrophysical systems like stars or inspiraling black holes, it is physically reason- 
able to assume that spacetime becomes flat in the limit of large distance, and that 
the systems can be considered "isolated", i.e. unaffected by the large scale struc- 
ture of the universe. The formulation of a rigorous concept of "asymptotic flatness" 
in GR is far from straightforward, due to the absence of a background metric or 
preferred coordinate system, in terms of which falloff rates can be specified. There 
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are two possible routes to overcome this problem: In the first approach one sim- 
ply assumes the existence of a suitable coordinate system, which is then used to 
formulate falloff conditions for tensor components in this coordinate system. We 
will follow this approach below to define quantities such as the total energy of a 
gravitating system. The obvious drawback is that taking limits is often problematic, 
in particular in a numerical context, and coordinate invariance of expressions has 
to be carefully checked. A resolution of these problems is provided by a definition 
of asymptotic flatness, where, after a suitable conformal rescaling of the metric, 
"points at infinity" are added to the manifold, one thus works on a compactified 
auxiliary manifold, and local differential geometry can be used to study the asymp- 
totic properties of the gravitational field [54]. Note that the notion of asymptotic 
flatness at timelike infinity does not make much sense in a general situation, because 
then all energy would have to be radiated away, leaving only flat space behind - 
excluding black holes or "stars" . The important notions are asymptotic flatness in 
spacelike and null (lightlike) directions. A detailed discussion of the global structure 
of spacetimes describing isolated systems in general relativity is not really possible 
here - an excellent overview is given in the textbook [23], for work in the context 
of numerical relativity see e.g. [55I56I57I58I27] , 

The concept of asymptotic flatness of isolated systems is intimately related to 
the possibility of defining the total energy-momentum for such systems in general 
relativity. In GR there exists no known well-defined local energy density of the 
gravitational field, but a total energy-momentum, which transforms as a 4-vector 
under asymptotic Lorentz transformations, can be assigned to an isolated system 
[23j , analogous to a particle in special relativity. It is a constant of motion and can 
therefore be expressed in terms of the initial data on an asymptotically flat Cauchy 
hypersurface. 

If a manifold has more than one asymptotically flat end, e.g. in the presence of 
wormholcs of the Einstcin-Rosen-bridge type, then - in general different - masses 
can be associated with each of these asymptotic regions. 

The expression for the energy momentum of GR at spatial infinity has been given 
first by Arnowitt, Deser and Misner in 1962 [29] in the context of the Hamiltonian 
formalism, and is usually called the ADM momentum, the time component being 
called the ADM energy or, somewhat inconsistently, the ADM mass, different from 
the rest mass to be defined below. 

The expressions for the mass and momentum are given as limits of surface 
integrals over non-covariant quantities, and have to be evaluated in asymptotically 
Cartesian (regular) coordinates {x 1 } - where the components of the metric tend to 
diag(l, 1,1) for large radii r — \J xix 1 . The surfaces are spheres S r of radius r. 

We define the surface integrals (which we will also refer to as ADM integrals) 

E(r) = ^J s V99 lJ 9 U (9ik,j ~ 9ij,k) dS h (27) 

p *(r) ^hjg^ 9 & ~ W dSt ' (28) 

h « = ^ji m J Vgx 1 (KL - K8l) dS t (29) 

which have to be evaluated in an asymptotically Cartesian coordinate system. 

The ADM energy Madm and linear and angular momentum Pj and Jj are then 
given by [59)28] 

M ADM = Urn E(r), Pj = lim P,-(r), Jj = lim J,(r), (30) 

r — ► oc r— >oo r — >oo 

and the rest mass Mr can be defined as M\ — M\ DM — Ylj=i 3 PjPj- 
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For an asymptotically Schwarzschildian metric h a b (here the extrinsic curva- 
ture falls off faster than in the general case , which allows for a boost, i.e. linear 
momentum) the ADM mass can be read off directly as the i-term of the metric: 

/ rn\ 4 const. , . 

hab ={ 1 + ^) °ab + Qab, | gab | < 1 + r 2 - ( 31 ) 

A fundamental issue of GR is the positivity of the ADM energy, since if the 
energy of an isolated system can be negative, it would most likely be unstable and 
decay to lower and lower energies. While it is trivial to write down a metric with 
negative mass, for reasonable matter fields with nonnegative energy density (thus 
satisfying the dominant energy condition), nonnegativity of the ADM energy thus 
can be expected on physical grounds. Indeed a complete proof of this positive energy 
conjecture has been given in 1982 by Schoen & Yau [T] (several simplified proofs 
have been given afterwards). 

For radiation processes we also require definitions of total energy, linear and an- 
gular momentum that decrease as energy and linear as well as angular momentum 
are radiated to infinity. The appropriate quantities are the Bondi quantities [60j . 
which can be defined as taking the limit of the ADM integrals not toward spatial 
infinity, but rather toward null infinity 61 62 63] . i.e., the limit to infinite distance 
is taken for constant retarded time instead of on a fixed Cauchy slice. In the con- 
text of our numerical treatment, the ADM and Bondi quantities can be calculated 
by computing values at several radii, and then performing a Richardson extrapola- 
tion (in extraction radius, not, as is more usual, in grid spacing). Here the Bondi 
quantities can be computed at any time for a fixed extraction radius, and have to 
be compared between different radii by taking into account the light travel time 
between the timelike cylinders of different radii, see e.g. [43] . 

Radiation quantities are conveniently defined in terms of the (complex) Bondi 
news function J\f(t) := dt{h + — Ih x ), which is the time derivative of the complex 
strain h taken at null infinity, where h+ and h x are the two polarization modes 
of the gravitational field, see eq. (Q}. The expressions for the radiated energy and 
momenta then become 



dE _ 1 

~dt ~ 167T 



NTJdn, (32) 

in 

T^zl Zi-NTJdQ, (33) 
ft 



dt 16tt 
dJ z 1 



where 



Re 

dt 16tt 



(- 



Ms / 77df dn 



(34) 



and an overbar denotes complex conjugation. The strain h is most often computed 
indirectly via double time integration of certain projections of the curvature tensor, 
see e.g. [43] for a detailed recent description in the context of numerical relativity. 



6 Horizons 



This section introduces the concept of trapped surfaces and apparent horizons. 
These are intimately related to two of the most fascinating features of general 
relativity: the appearance of singularities and the appearance of causal membranes: 
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event horizons that enclose so-called black holes - regions of spacetime that do not 
allow any information to escape to the outside world. 

For quite some time it was not clear whether singularities, which have first been 
found in highly symmetrical spacetimes such as the Friedmann-Robertson- Walker 
cosmologies or the Schwarzschild spacetime, are really generic, or mere artifacts 
of the high symmetry in these situations. Only in 1965 the generic character of 
singularities has been proven in a theorem by Penrose |64j . A variety of similar 
theorems with modified assumptions has followed, for an overview see e.g. Hawking 
and Ellis [22] or Wald [23]. 

A crucial notion in the singularity theorems is that of a trapped surface: a space- 
like 2-surface with the property that the area of the outgoing wavefront decreases 
toward the future. The central idea of the singularity theorems is then, that provided 
some energy condition {e.g. that the local energy density of matter is nonnegative) 
holds for the matter fields, gravity is always attractive, and light rays always get 
focused. Thus, if a region of spacetime starts to collapse, as is signaled by the ap- 
pearance of a trapped surface, this collapse cannot be halted and will continue until 
a spacetime singularity forms. 

The original Penrose singularity theorem [64] states that, provided a closed 
trapped surface exists and the Cauchy surface is non-compact, at least one of the 
following things will happen: 

1. There occurs negative local energy density. 

2. Einstein's equations are violated. 

3. The spacetime manifold is incomplete - a singularity occurs. 

4. The concept of spacetime loses its meaning at very high curvatures, possibly due 
to quantum gravity effects. 

This means, that in the framework of classical general relativity a trapped surface 
signals the occurrence of a singularity in the future. 

While the singularity theorems have shown that sufficiently "strong" initial data 
will develop a singularity, these theorems make no statement about the nature of 
the singularities. It is of particular interest, whether the singularities that arise in 
a physically reasonable collapse situation is visible to any observer. 

A wide-spread hope - at least in classical general relativity - is expressed by the 
cosmic-censorship hypothesis of Penrose |65) , which is by now floating through the 
literature in various formulations. The basic idea is that naked singularities - that 
is singularities that can be seen by outside observers - should not arise from regular 
initial data. The future light cone of the beginning of a naked singularity would be 
a Cauchy horizon and destroy predictability. A distinction is made between global 
and local nakedness. Globally naked singularities can influence asymptotic infinity 
(in an asymptotically flat setting) - they are not hidden by an event horizon. Inside 
of an event horizon, there may sit a locally naked singularity - while no information 
can escape from it to infinity it can still be seen by observers inside of the black hole. 
To rule out locally naked singularities all singularities have to be spacelike, which 
is equivalent to the spacetime being globally hyperbolic. In the weaker formulation, 
global hyperbolicity is only required outside of an event horizon. 

Counter-examples to overly restrictive formulations of cosmic-censorship are 
known in various cases, e.g. for phenomenological matter like dust outside of the 
physical validity of the equation of state, but also for sets of measure zero in the 
space of initial data (related to a collection of phenomena known as "critical col- 
lapse", see e.g. [66j for an extensive overview). 

Black holes, although dangerous and exotic, from this viewpoint are the good 
guys, the bad guys are naked singularities. From the viewpoint of quantum gravity, a 
violation of cosmic censorship actually seems rather desirable, since then regions of 
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strong curvature, where quantum effects are important, may actually be observable 
in the outside world. 

A notion that is derived from that of a trapped surface, is the apparent horizon: 
the boundary of the region of trapped surfaces, which turns out to be "marginally 
trapped" - i.e. the area of the outgoing wavefront neither increases nor decreases 
toward the future. In contrast to event horizons, which are defined in a global 
manner and can only be determined if the maximal time development is known, 
apparent horizons and trapped surfaces are defined locally in time, within a single 
slice, and can be directly determined from the initial data. Since the event horizon 
can only be defined if a suitable notion of infinity exists in a given spacetime, 
one could take the point of view of defining a black hole directly by means of 
apparent horizons 67J. The location in a given geometry and the analysis of the 
physical properties of the apparent horizon thus play an important role in numerical 
relativity. "Apparent horizon finders" are an essential component of all binary black 
hole evolution codes and discussed extensively in |68J. 



7 Initial data for numerical relativity 

7.1 The conformal approach to solve the constraints 

We have already taken a first step in identifying freely specifiable data for Einstein's 
equations by formulating them as a Cauchy problem. Given appropriate initial 
data on a spacelike hypersurface, the spacetime is determined as a unique time 
development (modulo diffcomorphisms) of these data for future and past times. Due 
to the presence of the constraint equations (fT2|) the initial data can not be specified 
freely, and the next task therefore is to extract the unconstrained part of the initial 
data in such a way that the constraints then uniquely determine the whole set of 
initial data, and thus the whole spacetime. There exists a standard formalism to 
accomplish this goal and solve the constraints, the conformal approach developed by 
Lichncrowicz, York, O Murchadha and others [59] • The conformal approach is based 
on conformal rescaling, in particular of the spatial metric h a b, which is expressed 
from a 'base metric' h a b and a conformal factor ip as 



By combination of ()35|) with a similar rescaling and decomposition of the extrinsic 
curvature into a traceless symmetric part and a part that is derived from a vec- 
tor potential, the constraints will be written as a coupled set of four well posed 
quasilinear elliptic PDEs for four 'gravitational potentials'. 

For the conformal rescaling (|35[) the scalar curvature transforms as 



h ab = ip 4 h a b, 



ip>0. 



(35) 



Next we define the trace-free part of the extrinsic curvature by 



= K 



ab 




The rescaling 



A ab = ip w A 



(36) 



then results in the property 



D a A ab = ^- w D a A a \ 



(37) 
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where D a is the unique derivative operator associated with the metric h a b- For K — 
const., eq. (|37p expresses the fact, that the momentum constraint is conformally 
invariant, if the tracefree part of the extrinsic curvature is properly rescaled. Note 
that the definition of A ab in eq. is equivalent to setting 

A ab = i> 2 A a b, 

where 

A a b = A cd h ac hbd, A a b = A cd h ac hbd- 
The traceless symmetric tensor A ab is decomposed as 

A ab = Af T + (LW) ab 

into a divergence-free (transverse) traceless part A^ T and a tracefree part that can 
be obtained from a potential W a , 

(LW) ab := D a W b + D b W a - -h ab D c W c . 

3 

Insertion of the reverse decomposition 

Af T = A ab - (LW) ab 
into the constraint equations (|12p yields 

D a {LW) ab = -D a A ab + \^D b K + 8Tnfj w j b , (38) 

- A s ^ + \Rh^ - l^ 7 (A ab - (LW) ab ) 2 + ^K 2 = (39) 

where we have now kept the energy density p and momentum density j b of matter 
fields in the expressions. The freely specifiable quantities here are the metric h ab , 
the trace of the extrinsic curvature, K, and a symmetric tracefree tensor A ab , which 
together comprise the local freedom in choosing initial data. The constraint equa- 
tions now take the form of a coupled elliptic system of PDEs for the 'potentials' ip 
and W a , the initial data are reconstructed as 

h ab = tfhab, K ab = (A ab - (LW) ab ) 2 %l)~ w + ^K^- A h ab . 

Prior to the choice of a set of fields (h a b, K, A ab ), one has to specify a 3- manifold 
S, on which the fields (h a b, K, A ab ) are defined and the equations (|39l38p are to be 
solved. If the manifold S has a nonempty boundary, or is not compact, it will be 
necessary to impose boundary conditions or asymptotic falloff conditions, which 
have to be chosen, along with the topology of S, on physical grounds. 

Due to the diffeomorphism invariance of GR, different initial data sets will give 
rise to the same spacetime (e.g. different slices of the same spacetime in different 
coordinate systems), which leads to the question for the number of local physical 
degrees of freedom represented by the initial data (h a b, K, A ab ). Since (h a b and K ab ) 
are both symmetric, the initial data are represented by 12 free functions. Three of 
these correspond to conditions on the spatial coordinate system (the metric h a b 
can be regarded as given by three functions by choosing a coordinate system where 
it is diagonal). In addition to initial data which are equivalent with respect to 
spatial diffeomorphisms, also data on different Cauchy surfaces may give rise to the 
same spacetime. A hypersurface is specified by one function, which can usefully be 
identified with the trace of the extrinsic curvature, also called the mean curvature - 
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since it corresponds to an average over the components of the extrinsic curvature. By 
imposing the 4 constraint equations on the remaining 8 free functions and dividing 
by two, one arrives at 2 local degrees of freedom per space point (confirming that 
GR is indeed a field theory). This is the same number of degrees of freedom as 
results from the linearized theory, which is the theory of a spin 2 field, where the 
2 degrees of freedom can be identified with 2 independent states of polarization - 
just as in the case of electromagnetism (see e.g. ref. [53], sec. 4.4). 

Most discussions of initial data restrict attention to so-called constant mean 
curvature (CMC) hypersurfaces. The condition that the mean curvature, or equiv- 
alently the trace of the extrinsic curvature, be constant is a coordinate independent 
statement and decouples the system (|39|38|) . The procedure of solution then be- 
comes the following: 

1. Choose h a bi A ab , K — const., 

2. solve dMJ for K ab , 

3. solve (|39|) for the scale factor ip, regarding the extrinsic curvature term as source. 

Slices with K = are called maximal slices. Maximal hypersurfaces embedded 
in a Lorentzian manifold locally maximize the 3-dimcnsional volume in the same 
way as a timelike geodesic maximizes the proper length, or a minimal 2-surface in 
Riemannian 3-space minimizes the area. The opposite extreme, i.e. a minimal 3- 
slice, or a maximal 2-surface is not possible, since deformations of the submanifold 
that just make the embedding 'more wiggly' will result in a change of volume of 
a definite sign: it is positive for a Riemannian hypersurface or a timelike curve in 
a Lorentzian manifold, and negative for a submanifold of a Riemannian manifold. 
The concept of extremal submanifolds thus is a natural generalization of geodesies 
as straightest curves. 

Considering the fact, that even simple metrics can be made to look arbitrarily 
complicated by coordinate choice, it is natural to specify initial data on hypersur- 
faces that are embedded as simply as possible, which leads to the consideration 
of maximal (ore more general CMC) slices as 'least wiggly'. Another advantage of 
maximal slices is that a foliation of maximal slices avoids singularities [70 , which 
has been utilized in many numerical calculations. 

Initial data for black holes can conveniently be constructed by "filling" the 
spacetime volume inside the horizon with artificial asymptotically flat ends. This 
construction enforces the presence of "throats" - minimal surfaces, and for non- 
vanishing extrinsic curvature also the presence of horizons (for vanishing extrinsic 
curvature the outermost minimal surfaces is actually an apparent horizon). These 
asymptotic regions are typically compactified for technical convenience, rendering 
the nontrivial topology of the resulting spacetime representable on R 3 , or S 3 if the 
physical asymptotically flat end is also compactified. Compactification naturally 
leads to singular behavior at the coordinate locations of the artificial asymptotic 
ends, which are commonly referred to as "punctures" . The treatment of the sin- 
gularity in the constraint equations is well understood [71] . The use of such initial 
data has first been advocated in [72], and has become the method of choice for a 
large fraction of work on binary black holes in numerical relativity following the 
prescription of Brandt and Briigmann [73] . The standard simplifying assumptions 
in the binary black hole literature are that the spatial geometry is conformally flat 
and that the extrinsic curvature is of the Bowen-York form [74], which is a family 
of solutions to the momentum constraint on a flat background, for which the total 
linear and angular momentum can be freely specified. By linearity of the momen- 
tum constraint and the flat background, superposition can be applied to construct 
multiple-black hole solutions. Removing the assumption of conformal flatness be- 
comes an issue particularly for spinning black holes, see [75] for a recent discussion. 
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An essential question regarding the construction of initial data is the issue of 
physical validity of the initial data set, say for the inspiral of two black holes. 
Data sets should contain little or no "artifical radiation", and should correspond 
to the actual inspiral of astrophysical compact objects. The typical aim is to start 
with initial data that correspond to the astrophysically most relevant case of quasi- 
circular inspiral, which essentially means that the orbital eccentricity is very small: 
eccentricity is radiated away rather efficiently for inspiraling black holes [7B]. It 
has recently been possible for the nonspinning equal mass case to directly use the 
initial momenta (and separations) computed from a PN inspiral, and show that not 
only the eccentricity of such initial data is very small, but also that the influence 
of "artificial radiation" inherent in the initial data can essentially be neglected. I 
expect this method to carry over to more general scenarios, for large spins of the 
black holes non-conformally flat initial data might have to be used along the lines 
of [75j and the references cited therein. For extensive reviews of the problem of 
construction initial data for general relativity see [77130) . 

8 The binary black hole revolution 

The inspiral and collision of compact objects - in particular of black hole binary 
systems - is considered one of the most important sources of gravitational waves for 
earth- and space-based detectors. Producing templates for gravitational-wave data 
analysis that describe signals from inspiraling compact binaries will require large 
parameter studies, and correspondingly large computational resources: The eventual 
goal of our simulations is to map the physical parameter space of gravitational wave 
signals from black hole coalescence, which is essentially given by the mass ratio and 
individual spins, as well as the initial orientation of the spins. The latter determines 
in particular the spin orientation at merger time, which may have a significant 
influence on the gravitational wave signal. 

In order to produce "complete" waveforms, which contain large numbers of grav- 
itational wave cycles from the inspiral phase, as well as the merger and ringdown 
phases, it is necessary to start the numerical simulations in the regime where Post- 
Newtonian analytical calculations are valid. These describe very accurately the 
waveforms of the early inspiral process, but break down for small separations of 
the black holes. This "matching" of analytical and numerical results requires large 
initial black-hole separations and large integration times. 

The numerical solution of the full Einstein equations represents a very complex 
problem, and for two black holes the spacetime singularities that are encountered in 
the interior of black holes pose an additional challenge. In order to obtain accurate 
results both the use of mesh refinement techniques and a good choice of coordinate 
gauge are essential. Together with the complicated structure of the equations - 
a typical code has between ten and several dozen evolution variables, and, when 
expanded, the right hand sides of the evolution equations have thousands of terms — 
this yields a computationally very complex and mathematically very subtle problem. 

For a long time, typical runs had been severely limited by the achievable evo- 
lution time before the simulations became too inaccurate or before the computer 
code became unstable, and there were serious doubts whether numerical relativity 
techniques could produce gravitational-wave templates, at least in the near future. 
This picture has drastically changed ever since in spring of 2005 Pretorius [78] 
presented the first simulation lasting for several orbits, using adaptive mesh refine- 
ment, second-order finite differencing, a sophisticated method to excise the singular 
interior of the black hole from the grid, and an implicit evolution algorithm. 

An alternative to the "excision" method of treating black holes is to "fill" the 
black hole with a topological defect in the form of an interior space-like asymp- 
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Fig. 1. Left panel: Coordinate tracks of the puncture location of one black hole for sim- 
ulations with grid sizes of the innermost boxes of {64 3 , 72 3 , 80 3 }, starting at a coordinate 
separation of D — 12M. Only in the last few orbits small differences are visible between 
the three runs discernable. Right panel: the waveform h+, rescaled with the extraction 
radius. Figure taken from |92] , 

totic end, the "puncture" [73], and to freeze the evolution of the asymptotic region 
through a judicious choice of coordinate gauge [79180140) . The latter approach, com- 
bined with a setup where the topological defect is allowed to move across the grid 
("moving puncture" approach |81|82J ) has led to a giant leap forward in the field 
[83I84I85I86I87I88I89| . taking the first orbit simulations of black holes [901 78j to 
more than ten orbits and allowing accurate wave extraction. 

In order to overcome phase inaccuracies in long evolutions, spectral methods 
have been suggested and significant progress has been made by the Caltech-Cornell 
group [91] . The Jena group has recently obtained good results with 6th order ac- 
curate finite differencing methods [32], which has enabled us to perform highly 
accurate evolutions over nine orbits (see fig. (J]) - units in the figures are given in 
units of the total initial black hole mass M) and perform a detailed comparison 
with predictions from PN approximations 92 44 18] • For an extensive overview on 
PN approximations see [T3] , 

A long standing problem has been the specification of non-eccentric initial data 
for numerical relativity simulations, which model inspiraling compact objects which 
have shed the eccentricity of their orbit through gravitational radiation. In [44j we 
show, that at least for the nonspinning equal mass case the initial momenta can be 
taken from a PN prescription. Part of the reason why this works is that the coor- 
dinate gauge used for numerical relativity agrees with the so-called ADMTT gauge 
adopted for PN calculations to excellent accuracy until very late in the inspiral, see 
fig. ([3|). The PN calculation is performed such that all initial eccentricity is shed 
during the first few hundred orbits, as shown in fig. ([5]). 

A particular focus of the last few months has been the so-called recoil or rocket 
effect due to "beamed" emission of gravitational radiation [9"4)95][96j . By momentum 
conservation, radiation of energy in a preferred direction corresponds to a loss of 
linear momentum and the black hole that results from the merger thus recoils from 
the center-of-mass frame with speeds of up to several thousand km/s. The velocity 
of this "kick" depends on the configuration of the system (e.g. the mass ratios and 
spins) and details of the merger dynamics, but not on the total mass (velocity is 
dimensionless in geometric units). 
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Fig. 2. The radial momentum component is plotted versus separation for PN-inspirals 
starting from D — 2QM and D = 40M. A separation of D = 20M is clearly not sufficient 
to produce non-eccentric inspiral parameters, since small oscillations can still be seen at 
D = 11M, while for D = 40M the initial eccentricity has essentially decayed away. Figure 
taken from |44] , 




Fig. 3. Orbital coordinate motion of a 9-orbits numerical relativity evolution compared 
with a PN evolution with the same initial parameters. In both panels the PN evolution is 
drawn as a dashed line. Top panel: the separation of the black holes (the puncture position 
in the full NR case). Bottom panel: the coordinate angular velocity. Figure taken from |18| . 



From an astrophysical point of view, the recoil effect is particularly interesting 
for massive black holes with masses > 10 5 M Q , which exist at the center of many 
galaxies and may have a substantial impact on the structure and formation of 
their host galaxies. The largest recoil effects have so far been found [93 97 for a 
particularly simple configuration: equal mass black holes with (initially) anti-aligned 
spins in the orbital plane. Such large kicks are on the order of 1% of the speed of 
light, and larger than the escape velocity of about 2000 km/s of giant elliptical 
galaxies. Smaller but still significant kick velocities have been found for several 
different types of black hole configurations [98|99|100|101|102|103j . 

Many challenges remain in the binary black hole problem: extreme mass ratios, 
combined, say, with large spins may remain problematic for some time; incorporat- 
ing realistic matter models adds a wealth of new problems. In order to perform large 
parameter studies, significant further optimizations for current codes are probably 
required, along with further mathematical insights and a better understanding of 
the general relativity aspects of the methodological foundations of the field. I be- 
lieve that with the "binary black hole revolution" the field of numerical relativity 
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Fig. 4. Coordinate positions of the black-hole punctures for model Mil from |93| up to 
t = 180. The black holes move out of the original plane and after merger the final black 
hole receives a kick in the negative z-direction. Figure taken from [93] . 

has found a new beginning rather than come to its end - and the study of colliding 
compact objects in particular will remain fruitful scientific ground for some time. 

I am grateful to Mark Hannam, Norbert Lages and Christof Gattringer for reading the 
manuscript and identifying some misprints and unclear points. 
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